NUMERICAL COMPUTATION OF SOLITON DYNAMICS 
FOR NLS EQUATIONS IN A DRIVING POTENTIAL 



MARCO CALIARI AND MARCO SQUASSINA 

Abstract. Wc provide some numerical computations for the soliton dynamics of the 
nonUnear Schrodinger equation with an external potential. After computing the ground 
state solution r of a related elliptic equation we show that, in the semi-classical regime, 
the center of mass of the solution with initial datum modelled on r is driven by the 
solution of X = —'W{x). Finally, wc provide some examples and analyze the numerical 
errors in the two dimensional case when V is an harmonic potential. 



1. Introduction 

1.1. Soliton dynamics behaviour. The goal of this paper is to provide a numerical 
investigation of the so-called soliton dynamics behaviour for the nonlinear Schrodinger 
equation with an external time independent smooth potential V 



iedt(j)e = -yA0, + Vix)(l), - X G M^, t > 0, 

0,(x,O) = 0o(x), xgM^, 



that is the qualitative behaviour of the solution ^^(t) in the semi-classical regime, namely 
for e (which of course plays the role of Planck's constant) going to zero, by taking as 
initial datum a (bump like) function of the form 

(/) M^) = r(^^)e^"-««, X G M^. 

We shall assume that A^>l,0<p<2/A^, iis the imaginary unit and r G fl C^(]R^) 
is the unique [Kwo] (up to translations) positive and radially symmetric solution of the 
elliptic problem 

(E) - ^Ar + \r = r'^^^^ in M^, 

for some value A > 0. Finally, Xq and are given vectors in that should be conveniently 
thought (in the transition from quantum to classical mechanics) as corresponding to the 
initial position and initial velocity respectively of a point particle. 
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In this framework, since (P) has a conservative nature, the typical expected behaviour 
is that the solution travels with the shape of r{{x — x{t))/e) (hence its support shrinks, as 
e gets small) along a suitable concentration line x{t) merely depending on the potential 
V and starting at xq with initial slope .^o- 

On the basis of the analytical results currently available in literature (see the discussion 
in Section 1.2), we believe that providing some numerical study is useful to complete 
the overall picture of this phenomenon and furnish some practical machinery for the 
computation of the solutions of {E) and, in turn, of {P)-{I). The authors are not aware 
of any other contribution in the literature on this issue. For the linear Schrodinger, some 
results can be found in [JY]. 



1.2. Facts from the theory. It is well-known that, given a positive real number m, the 
afore mentioned (ground state) solution r of [E] (where the value of A depends on m) 
can obtained through the following variational characterization on the sphere of L^(M^) 

(1.1) E{r) = mi{g{u) : u G H\R^), \\u\\l2 = m}, 
where S : H^{W^) — > M is the energy functional 

(1.2) S(u) = 1 [ Wul^dx — [ \u\^P+^dx. 



2 Jrn p 

Furthermore, there exists a suitable choice of m yielding A = 1 as eigenvalue in equa- 
tion (E). The restriction to the values of p below 2/A^ is strictly related to the global 
well-posedness of (P) for any choice of initial data 0o in H^. If p is larger than or 
equal to 2/N, then the solution can blow-up in finite time (see e.g. the monograph by T. 
Cazenave [Caz]). In particular, in the two dimensional case, p will be picked in (0, 1). 

From the analytical side, it has been rigorously known since 2000 that the solution 
(peit) of (P) remains close to the ground state r, in the sense stated here below, locally 
uniformly in time, as e is pushed to zero. As we said, this dynamical behaviour is typically 
known as soliton dynamics (for a recent general survey on solitons and their stability, see 
the work of T. Tao [Tao]). 

For the nonlinear equation (P), rigorous results about the soliton dynamics were ob- 
tained in various papers by J.C Bronski, R.L. Jerrard [BJ] and S. Keraani [Kee, Keel]. 
We also refer to [Squ] for a complete study of the problem with the additional presence 
of an external time independent magnetic vector potential A : M^, and to [MPS] 

for a study of a system of two coupled nonlinear Schrodinger equations, a topic which 
is rapidly spreading in the last few years. The arguments are mainly based upon the 
following ingredients: the energy convexity estimates proved by M. Weinstein [Wei, Weil] 
to get the so called modulational stability, the use of conservation laws (mass and energy) 
satisfied by the equation and by the associated Hamiltonian system in built upon the 
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guiding external potential V, that is the classical Newton law 

(x{t) = -VV{x{t)), 

(1.3) I x{0) = xo, 

U(o) = eo. 

Under reasonable assumptions on V (e.g. uniform boundedness of the second order partial 
derivatives), equation (1.3) admits a unique global solution (x(t),^(t)) which satisfies the 
following conservation law 

nt) = l\m\' + v{x{t)), n{t) = n{o), t>o. 

Let us now define a suitable scaling of the standard norm of H^{M.^) 

110111= ^'"^11 V0||i.+e-^||0||i., e>0. 

The precise statement of the soliton dynamics reads as follows 

Analytical Property 1.1 (Cf. e.g. [BJ, Keel]). Let (peit) be the solution to problem (P) 
corresponding to the initial datum (/) . Then there exists a family of shifts 6^ : 
[0, 27r) such that, as e goes to zero, 0e(x, t) is equal to the function 

(1.4) 0^(x,t) = X e M^, t > 0, 

up to an error function u!i;{x,t) such that \\i^e{t)\\m^ ^ ^(^)> locally uniformly in time. 

It is important to stress that, in the particular case of standing wave solutions of (P), 
namely special solutions of (P) of the form 

(j)e{x,t) =Us{x)e~^^\ xgM^, tGM+, {6 eR), 

where Us is a real-valued function, there is an enormous literature regarding the semi- 
classical limit for the corresponding elliptic equation 

See the recent book [AM] by A. Ambrosetti and A. Malchiodi and the references therein. 
To this regard notice that, if = (null initial velocity) and xq is a critical point of the 
potential V, as equation (1.3) admits the trivial solution x{t) = Xq and x{t) = for all 
t G M"*^, formula (1.4) reduces to 

0^(a;,t)=r(^^)e^''^W, x G M^, t > 0, 

so that the concentration of 0e(t) is static and takes place at xq, instead occurring along a 
smooth concentration curve in R^. This is consistent with the literature for the standing 
wave solutions mentioned above. 

For other achievements about the full dynamics of (P), see also [GSS, GSSl] (in the 
framework of orbital stability of standing waves) as well as [KN, KM] (in the framework 
of non-integrable perturbation of integrable systems). Similar results were investigated 
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in geometric optics by a different technique (WKB method), namely writing formally the 
solution as Ue = f/e(x, t)e'''*^^'*)/'^, with Ue = Uq + eUi + e^U2 ■ ■ ■ , where 9 and Uj are 
solutions, respectively, of a Hamilton- Jacobi type equation (the eikonal equation) and of 
a system of transport equations. In presence of a constant external potential, the orbital 
stability issue for problem (P) was investigated by T. Cazenave and P.L. Lions [CL], and 
by M. Weinstein in [Wei, Weil]. Then, A. Soffer and M. Weinstein proved in [SW] the 
asymptotic stability of nonlinear ground states of (P). See also the following important 
contributions: V. Buslaev and G. Perelman [BP], V. Buslaev and C. Sulem [BS], J. 
Frohhch, S. Gustafson, L. Jonsson, I.M. Sigal, T.-P. Tsai and H.-T. Yau [FGJS, FTY, 
ASPS], J. Holmer and Zworski [HZl], A. Soffer and M. Weinstein [SWl, SW2], T.-P. Tsai 
and H.-T. Yau [TY]. 

Another interesting problem concerns the case where the initial datum is multibump (for 
simplicity two bumps), say, 

(1.5) 0o(x) = ri(^^)ei^-«° +r2(^^)ei^-''o, x G M^. 
where are solutions to the problem 

for some fixed > 0, i = 1, 2 and xq, y^, C,o, f]o are taken as initial data for 

r x{t) = -VV{x{t)), ( y{t) = -VV{y{t)), 

(1.6) <^ a;(0) = xo, <^ y{0) = yo, 

Then we state the following 

Analytical Property 1.2 (Cf. e.g. [ASPS]). Let (j)£{t) he the solution to (P) correspond- 
ing to the initial datum (1.5). Then there exist two families of shifts 91 : 'K'^ [0,27r) 
such that, as e goes to zero, (j)s{x,t) is equal to the function 

(1.7) 0^(,;,t) =ri(^^^)eH-*W+^.^W] +r2(^^^)eM-^W^^^^ 

up to an error function ^^^(x, t) depending both on e and on the initial relative velocity 
V = \$^o — T]o\ (the larger is v the smaller is the error), locally uniformly in time. 

See figure 3 in the final section for a movie showing this behaviour. 



2. Numerical computation of the soliton dynamics 

In the numerical simulations included in the last section of the paper, we shall consider 
the two dimensional case. On the other hand, here we consider the general case. 
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^2 

iedt(l)e{x,t) = --A(l)e{x,t) + V{x)(l)e{x,t) - \(j)e{x,t)\^''(j)s{x,t) , X G M^^, 



2.1. Overview of the method. Our purpose is to solve the Schrodinger equation 
( 

ledtdeix, t) = 

(2.1) <^ ^ 2 

[ 0£(x, 0) = r£(x -xo), xe 
where r^{x) = u{x/e), so that (f>{x,t) = u{x)e~^^^ is the solution of 

(2.2) idt<i){x,t) = -^A0(x,t) - \(t){x,t)\''y{x,t) 

being u real, positive and minimizing the energy (1.2) under the constraint ||m||^2 
Instead of a direct minimization of the energy (see, e.g., [BT, CORT]), here we consider 
the following parabolic differential equation 

{dtr{x, t) = ^Ar(x, t) + r^P^\x, t) + A(r(x, t))r{x, t), x eR^, t>0 
r(x, 0) = ro(x), ||ro||^2 = m, x G 

with vanishing boundary conditions, where the map t i-h> A(r(-,t)) is defined by 

i 4^ I Vr(x, t) pdx - 4^ |r (x, t) p^+^dx 



m. 



A(r(x,t)) 



klli2 



This approach is similar to the imaginary time method (see, e.g., [BD]), based on the 
propagation of the Schrodinger equation along imaginary time —it and projection to the 
sphere of radius ^/rn. In equation (2.3), projection is not necessary and the energy 
decreases: in fact, if we multiply equation (2.3) by r(x, t) and integrate over M^, we easily 
get 

^) IIl2 = / r{x,t)dtr{x,t)dx = 



2 dt 

and if we multiply equation (2.3) by dtr{x,t) and integrate over M^, we get 

l^£{r{-,t)) = - [ \dtr{x,t)\^dx + \{r{x,t)) [ r{x,t)dtr{x,t)dx 
z at jj^N j^N 

\dtr{x,t)\'^dx < 



Hence, the steady-state solution roo(x) = r(x, t oo) of (2.3) satisfies ||roo||^2 = m 
and has a minimal energy. In fact, notice that by the results of [Kwo], for any A > 
there exists a unique (up to translations) positive and radially symmetric solution r = rx 
of (E). In turn, given Ai,A2 > 0, if ri,r2 : — > M denote, respectively, the positive 
radial solutions of the equations 

-^Ari + Airi = r^^+\ "^^^^^ + A2r2 = r^^^^, 
then it is readily verified that 



r2(x) = /iri(7x), 7 = ( ^ j , = f ^ 
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I 



- 0.05 



Figure 1 . The positive, radially symmetric and radially decreasing ground 
state solution Too of (1.1) with m = 1 and p = 0.2. Of course, in the compu- 
tation of Too , there is a spurious imaginary part of maximum value around 
10~^^, since the complex FFT algorithm is involved. The corresponding 
value of Aoo is Aoo = —0.37921. 



which tells us that that, up to a scaling, the solution corresponding to different values of 
A is unique. Notice now that, due to the choice of the bump like initial datum (Gaussian 
like, see (2.6)) in the iterations to compute Too (see the discussion below), it turns out 
that Aoo, defined as A(roo), is negative and r^o is positive, radially symmetric (see figure 1) 
and solves 

--Ar + A r = r^P+^ 

12 oo ~ '>oo' oo ' oo ' 

where Aoo = —Aoo > 0. If Tm denotes the ground state solution (with the corresponding 
positive eigenvalue denoted by A^), then we have 

(2.4) rm{x) = iirooi-fx), 7 = , ^ = (^^^ . 

On the other hand, by construction, we have 

namely /i^7~^ = 1. Finally, by the definition of 7 and // in (2.4), we get Xm = Aoo and 
7 = /i = 1, yielding from (2.4) the desired conclusion, that is 

''"00 ^m- 

Moreover, rao{x)e~^^^^°°^^^^^ is a solution of (2.2). We will take r^ix — xq) = roo{{x — 
XQ)/e) as our candidate initial condition for the time-dependent nonlinear Schrodinger 
equation (2.1). From a numerical point of view, it is convenient to compute directly 
fooix/e) instead of roo{x) and to apply the change of variable $(X, t) = v^e^0e(x, t), 
^/eX = X, to the nonlinear Schrodinger equation (2.1), and hence to equation (2.3). 
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Altogether, we need to solve 



(2.5) 
with 



dtR{X, t) = -AR{X, t) + e"^P/^R{X, tfP+^ + A(i?(X, t))R{X, t), X G 
R{X, 0) = Ro{X), WRoWh = me^, X e 

,J^^\VR{X,mX-e-^^/'J^^\R{X,t)\'^-^'dX 



A{R{X,t)) 



\R\\l2 



where R{X,t) = r{x/ e,t)\fe^ . Since it is not possible to numerically integrate the 
equation up to an infinite time, we will consider R{X, t) the steady-state as soon as 
S{R{X,t)) is stabilized within a prescribed tolerance. The initial condition Ro{X) can 
be arbitrarily chosen (in the class of bump like functions), but an initial solution with 
small energy will shorten the "steady-state" time i. Among the family of the Gaussian 
functions parameterized by a 



(2.6) R^iX) = V^a^/2e-l-^/v^lV2^(^^^ 
with lli^'^ll^a = me^ it is possible to choose the one with minimal energy. In fact 
S{R'') = a^- \VR\X)\MX -a^P- / \R\X)\^p+MX 

2 JrN P + 1 J^N 

If we define 

A=- \VR\X)\^dX, B = - / \R\X)\^^+^dX 

2 JrN P+l J^N 

the minimum for £{R") is attained for 



/ BNp\ 



The quantities A and B can be analytically computed and give 

N = l 

4 

N = 2 , B 



A 



2 ' 7r^P/2(p+ 1)1+^/2 

2"^^' N = 3 



4 

2.2. Numerical discretization. With the normalization introduced above, the nonlin- 
ear Schrodinger equation to solve is 

, idMX, t) = --A<^{X, t) + ^^^^^ $(X, t) - t), X G 

(2.7) ^ * ^ ' ^ 2 ^ ' ^ e ^ ^ y/^^ 

$(X,0) = i?(X-Xo,t), XgR^ 

A well-established numerical method for the cubic Schrodinger equation (focusing or de- 
focusing case) is the Strang splitting [BJM, BJM, CNT]. It is based on a split of the 
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full equation into two parts, in which the first is spectrally discretized in space and then 
exactly solved in time and the second has an analytical solution. We used the Strang 
splitting method as well. The first part is 

(2.8a) iai$i(X,t) = -iA$i(X,t) 

Thus, the Fourier coefficients of t) restricted to a sufficiently large space domain 

satisfy a linear and diagonal system of ODEs, which can be exactly solved. The second 
part is 

(2.8b) m.{X, t) = "^^^.(X, t) - ^±j^i,iX, t) 

It is easy to show that the quantity |$2(-^7^)P'' is constant in time for this equation. 
Then it has an analytical solution. Given the approximated solution ^ t„) 

of equation (2.7), a single time step of the Strang splitting Fourier spectral method can 
be summarized in 

(1) take $„(X) as initial solution at time t„ for (2.8a) and solve for a time step k/2, 
obtaining ^i{X,tn + k/2); 

(2) take <I>i(X, t„ + k/2) as initial solution at time tn for (2.8b) and solve for a time 
step k, obtaining $2(^5^n + k); 

(3) take $2(-^;^n + k) as initial solution for (2.8a) and solve for a time step k/2, 
obtaining $„+i(X). 

The result $„+i(X) is an approximation of $(X, t„ + k). Since the solutions of the first 
part and the second part are trivial to compute in the spectral space and in the real 
space, respectively, it is necessary to transform the solution from spectral space to real 
and from real space to spectral before and after step (2) above, respectively. All the 
transformations can be carried out by the FFT algorithm. The method turns out to be 
spectrally accurate in space and of the second order in time. 

Therefore, we used the Fourier spectral decomposition for the solution of equation (2.5), 
too. Together with the Galerkin method, it yields a nonlinear system of ODEs 



(2.9) 



R'{t) = '-DRit) + f{R{t)), t>0, 
^(0) = ^0, 



where R is the vector of Fourier coefficients, D the diagonal matrix of the eigenvalues of 
the Laplace operator and / the truncated Fourier expansion of the whole nonlinear part 
of equation (2.5). For the solution of equation (2.9) we used an exponential Runge-Kutta 
method of order two (see, e.g., [HO]), with the embedded exponential Euler method. 
Given the approximation /2„ ^ R{tn)i a single time step of the method is 

(1) set An+i = kn+i^D and Rni = -R„; 

(2) compute i?„2 = Gxp{An+i)Rni + kn+i'fi{An+i)f{Rni) (exponential Euler method); 

(3) compute Rn+l = Rn2 + kn+l'f2{An+l){-f{Rra) + f {Rnl)) 
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where (pi{z) and (p2{z) are the analytic functions 

6^ — 1 6^ — 1 — Z 

ipi{z) = , z ^0, ip2{z) = 5 , 2; 7^ 0, 

z z^ 

<^i(0) = 1, <^2(0) = \. 

The result is an approximation of R{tn + ^n+i)- Exponential integrators are explicit and 
do not suffer of time step restrictions. However, they require the computation of matrix 
functions. In our case, the matrices involved are diagonal and the computation of 

the matrix functions exp(/l„+i), <y9i(y4„+i) and (/)2(^n+i) is trivial. In order to compute 
the terms and f{Rn2), it is necessary to recover the functions in the real space 

corresponding to the Fourier spectral coefficients i?„i and Rn2, respectively, then compute 
the nonlinear part of equation (2.5) and finally to compute its Fourier transform. All the 
transformations can be carried out by the FFT algorithm. The term Rn+i—Rn2 in step (3) 
above can be used as an error estimate for R{tn+i) —Rn+i and then it is possible to derive 
a variable time step integrator. This is particularly useful for our aim of computing the 
steady-state of the equation: in fact, we expect that the as soon as the solution approaches 
the steady-state it is possible to enlarge the time step, thus reducing the computational 
cost. The method turns out to be spectrally accurate in space and of the second order in 
time. 



3. Two DIMENSIONAL EXAMPLES AND ERROR ANALYSIS 

In this section, in order to provide some examples, we reduce to the two dimensional 
setting and focus on the physically relevant case of harmonic potential 

V{x,y) = uJiX + LU2y , ci;i,ci;2>0, 

well-established in the theory of Bose-Einstein condensates. In the two movies starting 
from figure 2 we show the dynamics of the solitary wave along two Lissajous curves, 
periodic in the left side and ergodic for the right side. In the movie starting from figure 3 
we report the soliton dynamics in the case of an initial datum exhibiting a double bump 
behaviour (with a sufficiently large distance between the centers) up to the collision time. 
It is important to stress that in these figures the paths have an analytic expression and are 
plotted before the dynamics starts. The movies will then show that the centers of mass of 
the solitons follow adherently these curves up to the final computation time. An analysis 
of the error (in the single bump case) arising when the modulus of the solution |0£(x, t)| 
is replaced by the modulus of the expression in the representation formula (1.4), namely 
r{{x — x(t))/e), is indicated in figure 4. As predicted by the analytical property 1.1, the 
error in the || ■ Hh^ is below the order 0{e). 
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-4 -3 -2 -1 1 2 3 4 -4 -3 -2 -1 1 2 3 4 



Figure 2. In both simulation movies we set e = 0.01, p = 0.02, m = 1, 
{xo,yo) = (—3.0,-3.0), vq = (0,0). In the left movie, we chose ui = 1.4 
and ti;2 = 1 (rational ratio). In the right movie, we chose toi = and 
a;2 = 1 (irrational ratio). Notice that, although the ratios U2/0J1 are very 
close in the two examples, the soliton dynamics is ergodic in the right movie. 
Of course the figures refer to the (squared modulus of the) solution at the 
time t = and contain the concentration paths (admitting an analytic 
expression) that the soliton is going to travel on. 
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0.01 



0.001 



0.0001 



le-05 I ' ' ' ' ' ' — 

0.001 0.01 

e 

Figure 4. For the error analysis, we set p = 0.02, (xo,|/o) = (~0.5, —0.5), 
m = 1, u; = (2, 1) and a final time t = ir. With the change of variable 
we used, y/eX = x and U{X) = y/eu{x), we have ||^i||^2 = ||f^||£,2 and 
||Vxm||^2 = ^||Vx^||^2- Hence, the numerical error is computed through 
formula (written for the 2D case) \\u\\m^ = ^/^^^W^xUWj^i + ^~^||f^|li2- As 
predicted by the analytical property 1.1, the error in the || ■ Hh^ is below the 
order 0{e). 
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